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We report a search for gravitational waves from the inspiral, merger and ringdown of binary black 
holes (BBH) with total mass between 25 and 100 solar masses, in data taken at the LIGO and 
Virgo observatories between July 7, 2009 and October 20, 2010. The maximum sensitive distance 
of the detectors over this period for a (20,20)Mq coalescence was 300 Mpc. No gravitational wave 
signals were found. We thus report upper limits on the astrophysical coalescence rates of BBH as a 
function of the component masses for non-spinning components, and also evaluate the dependence 
of the search sensitivity on component spins aligned with the orbital angular momentum. We find 
an upper limit at 90% confidence on the coalescence rate of BBH with non-spinning components of 
mass between 19 and 28 Mq of 3.3 x 10~ 7 mergers Mpc _3 yr _1 . 
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I. OVERVIEW 

Binary black hole (BBH) systems are a major class of 
possible gravitational-wave (GW) sources accessible to 
ground-based interferometric detectors such as LIGO [1 
and Virgo [2J. As described in [3], for higher-mass BBH 
systems, the merger and ringdown stages of the coales- 
cence come into the detectors' sensitive frequency range 
and the search sensitivity is improved by using inspiral- 
merger-ringdown (IMR) matched filter templates. Such 
templates were used in [3J to search for compact binary 
coalescences (CBCs) signals with total masses between 
25 and 100 M Q in LIGO data. 

Our knowledge of possible high-mass BBH source sys- 
tems [1] is based on a combination of observations and 
astrophysical modelling: a summary of the recent evi- 
dence on both fronts is provided in [3j. A number of in- 
dicators point to the possibility of forming binary black 
holes with component masses mi , m 2 of ~ 20 — 30 M 
and beyond: in particular, predictions of the future fate 
of the high-mass Wolf-Rayet X-ray binaries IC10 X- 
1 and NGC 300 X-l 0; analyses of dynamical BBH 
formation in dense stellar environments; the growing 
evidence for the existence of intermediate-mass black 
holes (e.g. [6]); and population-synthesis modeling of 
low-metallicity environments [7 . A recent population- 
synthesis study [5] that considered a wide range of as- 
trophysical models found that in low-metallicity envi- 
ronments or under the assumption of weak wind-driven 
mass loss rates, the distribution of BBH chirp masses 
M = (mim 2 ) 3 / 5 (mi + m 2 ) -1 / 5 extended above 30 M Q ; 
for comparison, the chirp mass of a binary with compo- 
nent masses (50, 50) M is 43.5 M . 

In this paper we report a search for GW signals from 
coalescence of binary black holes with non-spinning com- 
ponents having masses mi, m 2 between 1 and 99 Mq 
and total mass M = mi + m 2 between 25 and 100 M Q , 
over the most recently taken coincident data from the 
LIGO and Virgo observatories. A companion paper [§] 
describes a search for low-mass binary inspiral signals 
with 2 < M/Mq < 25 over these data, while a search of 
previous LIGO and Virgo data for BBH mergers signals 
with total mass 100-450 M© is reported in [TU] . 

The joint science run used in this work — S6 at LIGO 
and VSR2 and VSR3 at Virgo — was the most sensitive to 
date to signals from coalescing BBH; this search was also 
the first for high-mass BBH coalescences in Virgo data. 
We describe the detectors and the joint S6-VSR2/VSR3 
science run in Section [TTJ The search pipeline used here is 
similar to that of [3[ [TT] , with changes to the ranking of 
events to account for variability of the noise background 
over the parameter space of the search and between de- 
tectors. We give a brief overview of the pipeline and de- 
scribe changes relative to previous searches in Sect ion [Hi] 

The output of the analysis is a set of coincident events 
where a potentially significant signal was seen in two 
or more detectors with consistent coalescence times and 
mass parameters. Events occurring at times when the de- 



tectors' environmental or instrumental monitor channels 
indicated a problem likely to corrupt the data are vetoed: 
either removed from the search or placed in a separate 
category, depending on the severity of the problem. The 
significance of each remaining candidate event is mea- 
sured by its false alarm rate (FAR) , the expected rate of 
noise events with a detection statistic value (defined in 
Section III I at least as large as the candidate's. 



As in previous LIGO- Virgo searches, the distribution 
of noise events in non-Gaussian data is estimated by ap- 
plying unphysical time-shifts to data from different detec- 
tors. Events with low estimated FAR are subject to a de- 
tailed followup procedure to check the consistency of the 
detector outputs around the event time and determine 
whether environmental disturbances or detector malfunc- 
tion could have caused a spurious signal at that time. 
The search did not find any significant gravitational- wave 
candidate events; we describe the most significant events 
in Section ITVl 

We then evaluated the sensitivity of the search to co- 
alescing BBH at astrophysical distances by analyzing a 
large number of simulated signals ("injections") added to 
the detector data. These are used to estimate the sensi- 
tivity of the search in terms of the sensitive distance in 
Mpc within which we would be able to detect a signal, 
averaged over the observation time and over source sky 
location and orientation, with significance above that of 
the loudest event observed in the search. From this we 
set upper limits on the rate of such coalescences as a 
function of their component masses. For this purpose we 
used two recently-developed families of IMR waveforms. 
The improved EOBNRv2 family [12] was used to assess 
the sensitive range of the search for comparison with the 
previous high-mass BBH search and to set upper lim- 
its on astrophysical coalescence rates; the IMRPhenomB 
waveform family [13] was used to assess the sensitivity of 
the search to coalescences of BBH with spinning compo- 
nents, where the component spins are aligned with the 
orbital angular momentum and thus the system does not 
precess. We describe the injections performed and the 
resulting sensitivity distances and upper limits in Sec- 
tion [v] 

To conclude, we briefly discuss outstanding issues for 
high-mass BBH searches and prospects for the advanced 
detector era in Section [VU 



II. S6 AND VSR2/VSR3 OBSERVATIONS 

The US-based Laser Interferometer Gravitational-wave 
Observatory (LIGO) comprises two sites: Hanford, WA 
and Livingston, LA. The data used in this search were 
taken during LIGO's sixth science run (S6), which took 
place between July 7, 2009 and October 20, 2010. Dur- 
ing S6 each of these sites operated a single 4km laser 
interferometer, denoted as HI and LI respectively. The 
2km H2 instrument at the Hanford site which operated 
in earlier science runs was not operational in S6. Follow- 
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ing LIGO's fifth science run (S5) Q], several hardware 
changes were made to the LIGO detectors in order to 
install and test prototypes of Advanced LIGO [JJ] tech- 
nology. These changes included the installation of higher- 
powered lasers, and the implementation of a DC readout 
system that included a new output mode cleaner on an 
Advanced LIGO seismic isolation table [15]. In addition, 
the hydraulic seismic isolation systems were improved by 
fine-tuning their feed-forward paths. 

The Virgo detector (denoted VI) is a single, 3km laser 
interferometer located in Cascina, Italy. The data used in 
this search were taken from both Virgo's second science 
run (VSR2) [16], which ran from July 7, 2009 to Jan- 
uary 11, 2010, and third science run (VSR3), which ran 
from August 11, 2010 to October 20, 2010. In the period 
between the first Virgo science run (VSR1) and VSR2, 
several enhancements were made to the Virgo detector: a 
more powerful laser and a thermal compensation system 
were installed, and noise due to scattered light in the out- 
put beams was studied and mitigated. During early 2010, 
monolithic suspensions were installed, which involved re- 
placing Virgo's test masses with new mirrors hung from 
fused-silica fibers. VSR3 followed this upgrade. 

Typical detector noise spectral densities as a function 
of frequency for S6 and VSR2/3 were given in [9 . Due 
to the improved low-frequency sensitivity of Virgo in 
VSR2/3 [T7], the lower frequency cutoff of our analy- 
sis for VI data was reduced to 30 Hz, compared to 40 Hz 
for LIGO data. 



III. DATA ANALYSIS PIPELINE 

Our search algorithm, which was described in detail in 
[3 [IT], is based on matched filtering the data in each 
detector against a template bank of IMR waveforms, 
recording local maxima of signal-to-noise ratio (SNR) as 
triggers, then testing these triggers for consistency of co- 
alescence time and mass parameters between two or more 
detectors via a coincidence test [18] , and for their consis- 
tency with the template waveform via the x 2 test [TS]. 
The x 2 test i s necessary to suppress noise transients (see 
Section III C ) , which cause a much larger rate of triggers 



with high SNR than expected in Gaussian noise. 



A. Filter templates and optimal search sensitivity 

As filter templates we used the same family of wave- 
forms as described in [3] constructed using the results of 
[2D] , which we will refer to as EOBNRvl. The parame- 
ter space covered by our templates was also unchanged, 
ranging from 1 to 99 M for the binary component masses 
mi, 7712, and from 25 to 100 M© for the total binary mass 
M = mi + rri2\ a search covering total mass values be- 
tween 2 and 25 M was reported in [9]. 

The recently-implemented EOBNRv2 [TJ] and IMR- 
PhenomB 113 waveforms are more accurate than their 
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FIG. 1: Horizon distances for non-spinning equal-mass IMR 
signals in the LIGO and Virgo detectors, using EOBNRv2 as 
signal model, averaged over periods of data when the detec- 
tor sensitivities were near optimal for S6 and VSR2 and -3, 
respectively. 



predecessors in that they are better fits to waveforms 
produced by numerical relativity (NR) simulations, and 
because the NR waveforms themselves have improved in 
accuracy and cover a wider range of parameter space 
[21] [22], [23, 24 . We investigated whether these improved 
waveform families could be efficiently detected by a bank 
of EOBNRvl filter templates. The method used was 
to calculate the overlap between an EOBNRv2 or IMR- 
PhcnomB signal waveform and an EOBNRvl template, 
maximizing over the parameters of the EOBNRvl tem- 
plate (see [35J for the general method). For EOBNRv2 
signals, over the range of mass ratio 1 < q = ml/m2 < 6 
and total mass 25 < M/M© < 100, in the worst case 
the maximized overlap (effectualness) of EOBNRvl tem- 
plates was greater than 0.97. For non-spinning IMR- 
PhcnomB signals, the smallest effectualness was greater 
than 0.98. Thus, the use of EOBNRvl templates did 
not significantly degrade the efficiency of our search, for 
non-spinning signals. The more recent waveform families 
are, however, useful in more accurately determining the 
sensitivity of the search to astrophysical BBH mergers. 

As a simple measure of the maximum possible search 
sensitivity, we show in Figure [TJ the horizon distances for 
equal-mass EOBNRv2 signals, defined as the distances 
at which an optimally-oriented coalescence directly over- 
head from a given detector would have an expected SNR 
of 8 for an optimal matched filter; as in [5] we average 
these distances over periods of data for which the detec- 
tor sensitivities were optimal, or close to optimal, over 
each science run. We see that the maximum sensitive 
distance for a (20,20)M© coalescence in each of the LIGO 
detectors was approximately 300 Mpc. The detector sen- 
sitivities varied significantly over the observation time of 
this search, as detailed in [26] , Sensitive distances aver- 
aged over observation time and over source sky location 



and orientation are reported in Section V B 
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FIG. 2: Representative distributions of SNR and %r values for simulated signal (red circle) and background (black ' x ') triggers 
in the LIGO detectors, with contours of the detection statistics used in the search. Note the systematically lower values of 
X 2 for background events with SNR p > 10 in shorter-duration templates (right plot) compared to longer-duration (left plot). 
Left — triggers with template duration greater than 0.2 s; dashed lines indicate contours of constant re-weighted SNR statistic, 
Eq. (IT]). Right — triggers with template duration below 0.2s; dashed lines indicate contours of constant effective SNR, Eq. B. 



B. Background estimate and event ranking statistic 

After obtaining a list of candidate coincident events, 
each consisting of two or more triggers with consistent 
template masses and coalescence times, we estimate the 
significance of each event relative to background. Our 
background distribution is obtained by finding coincident 
events after applying unphysical time-shifts (greater than 
light travel time) to data from different detectors; we 
performed 100 time-shifted analyses using multiples of a 
(0, 5, 10) s time offset for (HI, LI, VI) data, respectively. 
In order to compare the coincident events to background 
we require a ranking statistic. 

The aim of our ranking statistic is to optimize the sep- 
aration of signal from background in the search. We 
tuned it by studying the distributions of triggers over 
SNR (p) and x 2 for time-shifted background events, and 
for events resulting from simulated IMR signals added 
to the data ("software injections"). We found that the 
distribution of background triggers depended strongly on 
template duration, which varied between approximately 
0.05 s for templates with the highest total mass 100 M Q , 
in the LIGO detectors, to several seconds for lower-mass 
signals. We compare signal-background separation for 
triggers from longer vs. shorter-duration templates, in a 
representative period of LIGO data, in Figure[2j The per- 
formance of the x 2 test was markedly worse for templates 
shorter than 0.2 s in LIGO and VSR3 data, with some 
noise triggers in these templates having large SNR but a 
relatively small x 2 value, comparable to simulated IMR 
signals. For signals seen in LIGO, this threshold value of 
0.2 s corresponds to a total mass of approximately 45 M 
for equal-mass systems, or a total mass of approximately 
90 Mq for the most asymmetric templates used in the 
search. 

The poor performance of the x 2 test for short tem- 
plates in LIGO and VSR3 data can be attributed to the 



small number of template cycles over which the signal 
SNR is accumulated. By contrast, the x 2 test was found 
to be effective in penalizing noise artefacts in VSR2 data 
over the entire template bank of the search. 

We divided coincident events into two bins: one for 
which all participating triggers from HI or LI (or VI, 
in VSR3 data) had template durations above 0.2s ("long 
duration events" ) and one where at least one trigger from 
HI or LI (or VI, in VSR3) had a template duration below 
0.2s ("short duration events"). 

For all triggers participating in long duration events, 
and all VI triggers in VSR2 data, we used the re-weighted 
SNR p statistic [5] defined as 



p=<M(i + (x 2 ) 3 )/2] 1/6 Ar ' (1) 

[p for x 2 < 1, 

where x 2 = X 2 /(2p — 2), and where we chose the number 
of frequency intervals p used in the evaluation of x 2 |19j 
to be 10 [3J. For HI or LI triggers, or VI triggers in VSR3 
data, participating in short duration events we used the 
effective SNR statistic of |3J: 

The detection statistic, "combined SNR" p c , is then given 
by the quadrature sum of single-detector statistics, over 
the coincident triggers participating in an event. 1 

We calculate the detection statistic values separately 
in different coincident times (times when two or more 



Note that, as in previous searches, wc applied a weak p-de- 
pendent cut on the Xr values, to remove triggers with very low 
statistic values: the effect of this cut may be seen in Figure [2] 
where the top left of each plot is empty. 



detectors are recording data, labelled by the active de- 
tectors), due to their different background event distri- 
butions and different sensitivity to astrophysical signals. 
The FARs of coincident events are estimated by first com- 
paring their p c values to those of time-shifted background 
events in the same bin by duration, and with the same 
event type (i.e. the same detectors participating in coin- 
cidence) . The final detection statistic of the search, com- 
bined FAR, is determined by ranking the event's FAR 
against the total distribution of background FAR values, 
summed over both bins in template duration and over all 
event types within each coincident time: see Eq. (III. 7-8) 
of 0. 



C. Data quality vetoes 

The gravitational-wave strain data from the detec- 
tors contains a larger number of transient noise events 
(glitches) with high amplitude than would occur in col- 
ored Gaussian noise. In order to diagnose and remove 
these transients, each of the LIGO and Virgo observa- 
tories is equipped with a system of environmental and 
instrumental monitors that have a negligible sensitiv- 
ity to gravitational waves but may be sensitive to glitch 
sources. These sensors were used to identify times when 
the detector output was potentially corrupted [27H3D]. 
We grouped these times into two categories: periods with 
strong and well-understood couplings between non-GW 
transient noise sources and detector output, and peri- 
ods when a statistical correlation was found although a 
coupling mechanism was not identified. In our primary 
search, both for the identification of GW candidates and 
the calculation of upper limits, times in both these cat- 
egories, and any coincident events falling in these times, 
were removed ("vetoed") from the analysis. We also per- 
formed a secondary search for possible loud candidate 
events, in which only times with clear coupling of non- 
GW transients to detector output were vetoed. The total 
time searched for GW candidate events, in which only the 
first category of vetoes were applied, was 0.53 yr. 

Even after applying vetoes based on auxiliary (envi- 
ronmental and instrumental) sensors, significant num- 
bers of delta-function-like glitches with large amplitude 
remained unvetoed in the LIGO detectors. It was found 
that these caused artifacts in the matched filter output 
over a short time surrounding the glitch: thus, 8 s of time 
on either side of any matched filter SNR exceeding 250 
was additionally vetoed. Times removed from the pri- 
mary search by this veto were still examined for loud 
candidate events. 

Approximately 0.47 yr of coincident search time re- 
mained after applying all vetoes. Additionally, approxi- 
mately 10% of the data, designated playground, was used 
for tuning and data quality investigations. These data 
were searched for gravitational waves, but not used in cal- 
culating upper limits. After all vetoes were applied and 
playground time excluded, there was 0.09 yr of H1L1V1 



io 3 




Inverse false alarm rate (yr) 



FIG. 3: Cumulative distribution of coincident events found in 
the search vs. estimated inverse false alarm rate (FAR), over 
the total time searched for possible GW candidates, 0.53 yr. 
Grey contour shading indicates the consistency at la (dark) 
through 5a (light) level of search coincident events with the 
expected background. 

coincident time, 0.17 yr of H1L1 time, 0.10 yr of H1V1 
time, and 0.07 yr of L1V1 time, giving a total analysis 
time of 0.42 yr. 



IV. SEARCH RESULTS 

We found no significant or plausible gravitational-wave 
detection candidates above the noise background of the 
search. The cumulative distribution of coincident events 
found in the search vs. estimated inverse false alarm rate 
(FAR) is shown in Figure [3j The distribution is consis- 
tent with the expected background over the total time 
searched for GW candidates, 0.53 yr. 

The most significant coincident event found in the 
search, with lowest estimated FAR (highest inverse 
FAR), was at GPS time 939789782 and had an esti- 
mated FAR of 0.41 yr -1 . This event, an H1V1 coin- 
cidence in H1V1 coincident time with SNR values of 
31 in HI and 5.50 in VI and a combined SNR statis- 
tic value of p c = 11.98, occurred at a time when sev- 
eral short-duration, non-Gaussian transient excess power 
events were visible, over a period of several seconds, in 
time- frequency maps of the GW strain channel in HI. 
The high trigger SNR in HI was caused by the first of 
these transients; however the behaviour of the strain and 
SNR time series over the following few seconds is strongly 
inconsistent with a high-mass binary coalescence signal. 
The x 2 test value in HI was high (approximately 190), 
but not sufficiently high to rule out the trigger as a can- 
didate. 

The frequency spectrum of the noise transients in HI 
indicates their probable origin as stray light scattered 
into the interferometer beam. The appearance of the 
event in VI was consistent with both a quiet GW signal 
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FIG. 4: Cumulative distributions of coincident events and estimated background over combined SNR statistic p c , over the total 
time searched for possible GW candidates. Grey shaded bands indicate la— 5a consistency with the estimated background 
distribution. Left — Distribution for "long duration" events. Right — Distribution for "short duration" events. The two event 
bins and the combined SNR statistic are described in Section [Til 13 



and with random noise. The event time was not vetoed in 
the search, as none of the instrumental or environmental 
channels at HI which we found to be significantly corre- 
lated with noise transients in the GW channel showed a 
malfunction at the event time. 

The second loudest event was at GPS 963363786 with 
an estimated FAR of l.Oyr -1 . It was an H1L1 coinci- 
dence in H1L1 time with SNRs of 13 in HI and 70 in 
LI and a combined SNR statistic value of p c = 10.48. 
This time was subject to a veto in HI, due to a problem 
in a high-voltage power supply near the output photodi- 
odes (affecting some weeks of data) which caused bursts 
of broad-band non-Gaussian noise. However, since this 
excess noise was not sufficiently severe to preclude detec- 
tion, the time was still searched for possible high-SNR 
candidate events, as described in Section |ill C| The trig- 
ger in LI was caused by a high amplitude non-Gaussian 
transient of very short duration, part of a population 
of sporadic glitches for which no effective veto could be 
found. This event also failed a detailed followup, as time- 
frequency maps of excess power, and the time series of 
SNR and \ 2 in HI, were inconsistent with an IMR signal. 

The next few most significant events had estimated 
FARs of a few per year and were thus entirely consistent 
with background. 

We show the p c distributions of coincident events 
and time-shifted background events, for the total time 
searched for possible GW candidates, in Figure |4j sepa- 
rating the long- and short-duration events since p c is a 
different function of p and \ 2 m each. The p c values of all 
the loudest search events were less than 12 (for compari- 
son, a BBH coalescence signal with an SNR of 8 in each 
of two detectors would give approximately p c ~ 11.3). 
Thus, the data quality veto procedure, in combination 
with the signal-based x 2 test, were sufficient to remove 
or suppress loud detector artefacts to a level where the 
sensitivity of our search is not greatly impaired. 



V. UPPER LIMITS ON BBH COALESCENCE 
RATES 

Given the null result of our search for BBH coalescence 
signals, we wish to set observational limits on the astro- 
physical rates of such signals, over regions of parame- 
ter space where our search has non-negligible sensitivity. 
As discussed in [5] , the distance reach of this high-mass 
search is such that the source population can be treated 
as approximately homogeneous over spatial volume; thus 
we aim to set limits on the rate density of coalescences, 
in units of Mpc _3 yr _1 . 



A. Upper limit calculation procedure 

Our upper limit calculation is similar to that per- 
formed in [9] : it is based on the loudest event statistic [31] 
applied as described in Section V of [3] with minor im- 
provements in implementation. We divide the data into 
9 periods of approximately 6 weeks each; in 4 of these 
only H1L1 coincident time was recorded, whereas in the 
remaining 5 we had four types of coincident time (H1L1, 
H1V1, L1V1 and H1L1V1). 

In each of the resulting 24 analysis times, we estimate 
the volume to which the search is sensitive by reanalyz- 
ing the data with the addition of a large number of simu- 
lated signals ( "software injections" ) in order to model the 
source population. Our ability to detect a signal depends 
upon the parameters of the source, including the compo- 
nent masses and spins (magnitudes and directions), the 
distance to the binary, its sky location, and its orien- 
tation with respect to the detectors. Numerous signals 
with randomly chosen parameter values were therefore 
injected into the data. 

To compute the sensitive volume over a given range of 
binary masses ("mass bin"), we perform a Monte Carlo 
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integration over the other parameters to obtain the effi- 
ciency of the search — determined by the fraction of sim- 
ulated signals found louder than the loudest observed 
coincident event in each analysis time — as a function of 
distance. Integrating the efficiency over distance then 
gives the sensitive volume for that analysis time, and an 
associated sensitive distance. 

We then estimate the likelihood parameter A of signal 
vs. background at the combined FAR value of the loud- 
est observed event in each analysis time, for each mass 
bin, as described in [3l (3Tj . Using these A values and 
the estimated sensitive volumes we find the probability 
of the measured loudest FAR value as a function of the 
astrophysical merger rate, i.e. the likelihood of the data 
given the model, in each analysis time and for each mass 
bin. Given a prior probability distribution over the rate, 
in each mass bin, we then multiply by the likelihoods of 
the loudest events from all the analysis times to form a 
posterior over rate: see Section V of [3], and [32], for 
relevant formulae. 

The likelihood function for each analysis time depends 
on the sensitive volumextime searched; however, the sen- 
sitive volumes have statistical uncertainties due to the 
finite number of injections performed in each mass bin, 
and systematic uncertainties in the amplitude calibra- 
tion for each detector. As detailed in [5], we take an 
overall 42% uncertainty in volume due to calibration er- 
rors. We marginalize over statistical uncertainties for 
each analysis time separately, but since systematic cal- 
ibration errors may be significantly correlated between 
analysis times we perform this marginalization |33j once 
after combining the likelihoods from all analysis times. 
We then find the 90% confidence upper limit based on 
the marginalized posterior distributions over rate. 

We also calculated an average sensitive distance in each 
mass bin, defined as the radius of a sphere such that 
the sphere's volume, multiplied by the total search time, 
equals the total sensitive volumextime over all analysis 
times. 

Since the injected waveforms are phenomenological 
models, our upper limits will also be systematically af- 
fected to the extent that the true IMR waveforms differ 
from these models. These uncertainties are difficult to 
quantify over the search as a whole and we will not at- 
tempt to incorporate them into our quoted limits. A 
comparison of EOBNRv2 waveforms against numerical 
relativity simulations for mass ratios q = 1,4,6 shows 
possible SNR biases of at most a few percent within the 
total mass range 25 < M/M < 100. 

The rate priors that we use for non-spinning EOB- 
NRv2 signals in the S6-VSR2/3 search are derived from 
the results of the S5 high-mass BBH search [3]. How- 
ever, the prior results from the S5 search were affected 
by an incorrect treatment of marginalization over errors 
in the sensitive volume and flaws in the numerical proce- 
dure used to estimate the A values, resulting overall in an 
over-conservative set of 90% rate upper limits [32] . These 
problems were recently addressed, leading to revised up- 



per limits from S5 data; the resulting revised posteriors 
over coalescence rate were used as priors for our main 
upper limit calculation. Revised rate upper limits from 
S5 data alone are included in Table|T] Note that, as priors 
for the S5 calculation, uniform probability distributions 
over rate was taken; this uniform prior is a conservative 
choice for setting upper limits. 



B. Rate limits from EOBNRv2 injections 

In order to evaluate the search sensitivity to non- 
spinning IMR signals over a wide parameter space, and 
to allow a comparison to previous search results, we used 
an implementation of the recently developed EOBNRv2 
waveform family [T2] as simulated signals. 

The EOBNRv2 waveform family was designed using 
results from [3~4Tf3"T] . The inspiral waveforms in the EOB- 
NRv2 model are improved over EOBNRvl by calibrating 
two adjustable parameters against five highly accurate 
NR simulations of mass ratios q = 1,2,3,4 and 6, gen- 
erated by the pseudospectral code SpEC [53J [21] . These 
two parameters are the pseudo-4PN and 5PN coefficients 
a§ and a§ entering the EOB radial potential (see, e.g., 
Eq. (II7a) in Ref. [3] and related discussion.) EOBNRv2 
also improves over EOBNRvl by including all known 
post-Newtonian (PN) corrections in the amplitude, by 
using a more accurate estimate of the radiated energy 
flux, by dropping the assumption of quasi-circular orbits, 
by improving the matching of the inspiral-plunge wave- 
form to the ringdown modes, and by improving the ex- 
trapolation to large mass ratios. The differences between 
these EOBNRv2 and NR waveforms are comparable with 
numerical errors in the NR simulations. 

EOBNRv2 injections were distributed to "over-cover" 
the parameter range of the search, in order to ensure 
complete coverage of the mass bins used in [5] . The in- 
jections were distributed approximately uniformly over 
the component masses mi and mi, within the ranges 
lM <m,< 99 M and 20 M Q < M < 109 M Q . 

The resulting 90% confidence upper limits on non- 
spinning coalescence rates are displayed in Figure [5] left 
plot, and in Table [i] These upper limits supersede those 
reported in ,3'. As explained in the previous section, we 
used revised priors over rate from the S5 analysis. These 
were obtained using (non-spinning) EOBNRvl simulated 
signals, and a smaller number of IMRPhenomA simula- 
tions. As described in [3] their distances were appro- 
priately adjusted to account for the amplitude bias of 
the older waveform families, by comparison with current 
NR simulations of BBH merger. Due to the restricted 
range of parameters of NR simulations available at the 
time of the previous analysis, we quote limits on astro- 
physical rates only for bins within the range of mass 
ratio 1 < q < 4. For binaries with both component 
masses lying between 19 and 28 M we find a 90% limit 
of 3.3 x 10~ 7 mergers Mpc~ 3 yr -1 . 

The averaged sensitive distance to simulated EOB- 
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FIG. 5: Left — Upper limits (90% confidence) on BBH coalescence rates in units of 10 _7 Mpc _3 yr _1 as a function of binary 
component masses, evaluated using EOBNRv2 waveforms. Right — Average sensitive distance for this search to binary systems 
described by EOBNRv2 signal waveforms, in Mpc. 



TABLE I: Search sensitive distances and coalescence rate upper limits, quoted over 9 M0-wide component mass bins labelled 
by their central values. We also quote the chirp mass M at the centre of each bin. The sensitive distance in Mpc (averaged 
over the observation time and over source sky location and orientation) are given for EOBNR waveforms in S5 data rescaled 
for consistency with NR results [3], and for EOBNRv2, IMRPhenomB non-spinning ("PhenomB nonspin") and IMRPhenomB 
spinning ("PhenomB spin") waveforms in the S6-VSR2/3 data. The last two columns report 90%-confidence rate upper limits 
in units of 10 -7 Mpc _3 yr _1 , for bins with component mass ratios 1 < mi/7712 < 4, for S5 data (revised relative to j3j) and the 
cumulative upper limits over S5 and S6-VSR2/3 data, as presented in this work. 
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Merger Rate (Mpc~ 3 yr _1 ) 

FIG. 6: Cumulative posterior probabilities over astrophysical 
merger rate, for the bins shown in Figure[5]with central values 
mi = m 2 = 50, 41, 32, 23, 14 M (left to right). We show 
the probability level corresponding to the 90% confidence rate 
limit (dashed horizontal line). These posteriors were evalu- 
ated for signals described by the EOBNRv2 waveform family 
in S6 data using S5 search results as prior information. 




Total mass (M Q ) 



FIG. 7: Dependence on aligned spin and total mass of the 
averaged sensitive distance of our search to phenomenological 
inspiral-merger-ringdown waveforms. For each of 6 bins in to- 
tal mass M, we show the sensitivity for IMRPhenomB signals 
with negative aligned spin parameter \ (left), non-spinning 
signals (centre) and signals with positive aligned spin param- 
eter (right). The simulated signal parameters were restricted 
to mass ratios 1 < q < 4 and aligned spins —0.85 < \ < 0.85. 



NRv2 waveforms over the S6-VSR2/3 observation time, 
for the entire parameter space of the search, is displayed 
in Figure [5j right plot. 

To illustrate our statistical upper limit method we dis- 
play the cumulative posterior probabilities over coales- 
cence rate, for a selection of the mass bins of Figure [5] 
covering the equal-mass line mj — rn 2 , in Figure [6] This 
figure shows the dependence of the quoted upper limit on 
the confidence level; we choose to use a 90% confidence 
limit as indicated by the dashed line. 



C. Sensitivity to non-spinning and spinning 
IMRPhenomB injections 

The spins of the component black holes are known to 
have a potentially large effect on the emitted IMR wave- 
form (e.g. [Ml EH]): affecting the phase evolution and, in 
the case of spins significantly out of alignment with the 
orbital angular momentum, producing amplitude modu- 
lations due to precession [50] . X-ray observations of the 
spins of accreting black holes in binary systems, while 
technically challenging, indicate a fairly uniform distri- 
bution over the entire range < a = S/m 2 < 1 [4TH47] , 
Note that such measurements apply to black holes in X- 
ray binaries, which may not be representative of spins in 
BBH systems. 

Indications that spin-orbit misalignment in field bina- 
ries may be small come from observations of the micro- 
quasar XTE J1550-564 [35], and from population synthe- 
sis models |49j . For dynamically formed binaries, how- 
ever, the component spins may be largely independent of 
each other and of the orbital parameters. 

In any case it is desirable to perform injections using 



spinning BBH coalescence waveforms, in order to esti- 
mate how far our search was sensitive to such signals. 
Knowledge of spinning BBH coalescence waveforms is, 
however, currently limited to a relatively small number 
of numerical relativity simulations (see [13. and refer- 
ences therein), most of which have spin magnitudes sig- 
nificantly below unity, and only a few including non- 
aligned spins. Only very recently have simulations for 
near-extreme spins [50] been performed. Thus, as a first 
step towards quantifying the sensitivity of the search to 
spinning waveforms over a broad parameter space, we use 
the IMRPhenomB waveform family 13] which models 
IMR signals from BBH with aligned/anti-aligned spins. 
This waveform family is parameterized by the total mass 
M = mi + mi-, the mass ratio q = mi/mi and a single 
aligned spin parameter %, 2 defined for spins Si parallel 
to the orbital angular momentum as 



mi 
M 



Xi 



m 2 
M 



X2, 



(3) 



where Xi = Sj • L/mf is the dimensionless spin of black 
hole i projected onto the orbital angular momentum 
L. These waveforms are calibrated against numerical- 
relativity simulations in the parameter range 1 < q < 4 
and —0.85 < x < 0.85, and, for the inspiral part, to 
the calculated evolution in the extreme-mass-ratio (test 
mass) limit. The waveform family is constructed in the 



2 This parameter is not to be confused with the x 2 test mentioned 
earlier in Section IIIII 
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frequency domain and then converted to the time domain 
by means of an inverse Fast Fourier Transform. 

As simulated signals we used two sets of IMRPhenomB 
injections, a non-spinning set and a spinning set. Both 
were uniformly distributed in total mass between 25 and 
100 Mq, and uniformly distributed in q/{q + 1) = mi/M 
for a given M, between the limits 1 < q < 4. In addi- 
tion, the spinning injections were assigned aligned spin 
components Xi uniformly distributed between -0.85 and 
0.85. 

To illustrate the effect of aligned spin on the search 
sensitivity, we plot in Figure [7] the average sensitive dis- 
tance over the S6-VSR2/3 observation time, in bins of 
total mass M, for both non-spinning simulated signals 
and for injections with \ < and x > respectively. 

Component spin is expected to have several effects on 
our search, compared to its performance for non-spinning 
systems. First, the amplitude of the expected signal 
from a coalescence at a given distance may depend on 
spin: see for instance Figure 3 of |13j . where the horizon 
distance for IMRPhenomB signals with optimally fitting 
templates, with Initial LIGO noise spectra, was found to 
increase steeply with increasing positive x- Second, the 
EOBNR templates used in our search may have reduced 
overlap with the simulated spinning signals, leading to a 
loss of sensitivity. Third, the signal-based x 2 test values 
are expected to be higher than if exactly matched spin- 
ning templates were used, due to "un-matched" excess 
power in the signals; this would further reduce the search 
sensitivity. Given the complexity of the search pipeline, 
it is not clear which effect would dominate. Figure [7J in- 
dicates higher sensitivity to positive- % signals even with 
the current non-spinning templates, but also shows that 
the search is significantly less sensitive to negative- x sig- 
nals at higher values of total mass M. 



D. Waveforms including higher spherical harmonic 
modes 

In the filter templates and in all injections used in this 
search, we consider only the dominant mode of GW emis- 
sion from coalescing binaries, the (l,m) — (2,2) mode. 
In general, higher-order modes are important in BBH 
with asymmetric masses and significant component spins. 
Omitting these modes in templates will neglect their con- 
tributions to the SNR [5TJ [55] , and may also lead to a 
worse (higher) value of the % 2 test, tending to reduce 
the sensitivity of the search. However, such effects will 
depend strongly on the mass ratio and on extrinsic (an- 
gular) parameters, and it is beyond the scope of this anal- 
ysis to investigate them in detail. In some quanti- 
tative calculations of mismatch between templates with 
and without higher modes were made, for a few sets of 
parameter values, using one projected Advanced LIGO 
sensitivity curve. 



E. Astrophysical implications 

There is no commonly accepted astrophysical rate es- 
timate for high-mass CBCs, owing to the many pos- 
sible formation scenarios and the considerable uncer- 
tainties affecting them. In [5J [S3] a rate of 3.6l2g x 
10~ 7 Mpc _3 yr _1 for IMR signals from binaries with chirp 
mass comparable to 15 M Q was estimated based on two 
observed tight binaries believed to consist of a mas- 
sive stellar BH and a Wolf-Rayet star. Our 90% up- 
per rate limit for the bin with component masses 19 < 
mi/M < 28, 10 < m 2 /M Q < 19, for which the chirp 
mass ranges between 12 and 20 M Q , is 5.9 x 10~ 7 merg- 
ers Mpc~ 3 yr -1 . Thus, current searches are close to the 
sensitivity necessary to put nontrivial constraints on as- 
trophysical scenarios of BBH formation and evolution. 
However, we remind the reader that systems with near- 
extremal (Si/mf > 0.85) or significantly non-aligned 
spins, or for which higher signal harmonics make a con- 
siderable contribution to the waveform seen at the detec- 
tors, were not included in our sensitivity studies. 

VI. DISCUSSION 

The present search is an advance over that reported in 
[3] in three main respects: the improved sensitivity of the 
LIGO and Virgo detectors over previous science runs; im- 
proved understanding of the search background by identi- 
fying the duration of the IMR templates as the dominant 
parameter controlling their output in non-Gaussian de- 
tector data; and the use of updated, more accurate signal 
models to assess search sensitivity, including models de- 
scribing component spins aligned to the orbital angular 
momentum. 

There are, however, many issues that remain to be ad- 
dressed in order for future data from advanced detectors 
[T4l [S"4"] to be best exploited in searching for high-mass 
CBC. Among these, the metric currently used for tem- 
plate bank placement and for testing mass coincidence 
between detectors is based on the inspiral portion of CBC 
waveforms only. The search may be improved by using 
a more accurate metric for IMR waveforms, and more 
radically by also including spin-aligned IMR templates 
as matched filters, which may significantly increase the 
sensitivity to spinning BBH. 

Separating signal from non-Gaussian noise events in 
short filter templates, where our signal-based x 2 tests 
are not effective, remains a difficult problem. Improved 
methods, including the use of amplitude consistency tests 
between different detectors and multivariate classifiers, 
are currently being investigated. 

Advanced detectors coming on line in the coming years 
will improve the sensitivity to GW relative to the first 
generation by a factor 10 over a broad frequency range, 
and will achieve good sensitivity down to a low-frequency 
limit of ~ 10 Hz. The volume of space over which future 
searches will be sensitive to IMR signals is therefore ex- 
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pected to increase by a factor 10 3 or more depending on 
the binary masses. We thus expect to extract significant 
information on BBH source populations over the param- 
eter space of future searches. 
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